#*************        丁香园-NHANES     *************
#*************        Analyze Code      *************
#*************   RCS 非线性拟合项讲解   *************
#*************                          *************

#### 0.准备好环境 ####
library(gtsummary)
library(survey)
library(haven)
library(rms) # RCS-Regression Modeling Strategies
# 运行一下读取 mort 数据的函数
NHSMortRead <- function(file.path){
  srvyin <- file.path   # full .DAT name here
  srvyout <- "data" # shorthand dataset name here
  
  # Example syntax:
  #srvyin <- paste("NHANES_1999_2000_MORT_2019_PUBLIC.dat")   
  #srvyout <- "NHANES_1999_2000"      
  
  
  # read in the fixed-width format ASCII file
  dsn <- read_fwf(file=srvyin,
                  col_types = "iiiiiiii",
                  fwf_cols(seqn = c(1,6),
                           eligstat = c(15,15),
                           mortstat = c(16,16),
                           ucod_leading = c(17,19),
                           diabetes = c(20,20),
                           hyperten = c(21,21),
                           permth_int = c(43,45),
                           permth_exm = c(46,48)
                  ),
                  na = c("", ".")
  )
  
  # NOTE:   SEQN is the unique ID for NHANES.
  
  # Structure and contents of data
  str(dsn)
  
  
  # Variable frequencies
  
  #ELIGSTAT: Eligibility Status for Mortality Follow-up
  table(dsn$eligstat)
  #1 = "Eligible"
  #2 = "Under age 18, not available for public release"
  #3 = "Ineligible"
  
  #MORTSTAT: Final Mortality Status
  table(dsn$mortstat, useNA="ifany")
  # 0 = Assumed alive
  # 1 = Assumed deceased
  # <NA> = Ineligible or under age 18
  
  #UCOD_LEADING: Underlying Cause of Death: Recode
  table(dsn$ucod_leading, useNA="ifany")
  # 1 = Diseases of heart (I00-I09, I11, I13, I20-I51)
  # 2 = Malignant neoplasms (C00-C97)
  # 3 = Chronic lower respiratory diseases (J40-J47)
  # 4 = Accidents (unintentional injuries) (V01-X59, Y85-Y86)
  # 5 = Cerebrovascular diseases (I60-I69)
  # 6 = Alzheimer's disease (G30)
  # 7 = Diabetes mellitus (E10-E14)
  # 8 = Influenza and pneumonia (J09-J18)
  # 9 = Nephritis, nephrotic syndrome and nephrosis (N00-N07, N17-N19, N25-N27)
  # 10 = All other causes (residual)
  # <NA> = Ineligible, under age 18, assumed alive, or no cause of death data available
  
  #DIABETES: Diabetes Flag from Multiple Cause of Death (MCOD)
  table(dsn$diabetes, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  #HYPERTEN: Hypertension Flag from Multiple Cause of Death (MCOD)
  table(dsn$hyperten, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  # Re-name the dataset, DSN, to the short survey name then remove other R objects
  assign(paste0(srvyout), dsn)
  rm(dsn, srvyin, srvyout)
  # View(data)
  return(data)
}


#### 1. RCS 建模（以 Logistics 为例） ####

##### 1.1 构建数据集 ####
n <- 1000 # define sample size set.seed(17) # so can reproduce the results
treat <- factor(sample(c('a','b','c'), n, TRUE))
num.diseases <- sample(0:4, n, TRUE) 
age <- rnorm(n, (50)^(1/2), 10)  # 后续设置 age 的平方项
cholesterol <- rnorm(n, 200, 25)
 weight <- rnorm(n, 150, 20) 
sex <- factor(sample(c('female','male'), n, TRUE)) 
label(age) <- 'Age' # label is in Hmisc 
label(num.diseases) <- 'Number of Comorbid Diseases' 
label(cholesterol) <- 'Total Cholesterol' 
label(weight) <- 'Weight, lbs.' 
label(sex) <- 'Sex' 
units(cholesterol) <- 'mg/dl' # uses units.default in Hmisc 
# Specify population model for log odds that Y=1
L <- .1*(num.diseases-2) + .045*(age^2-50) +
  (log(cholesterol - 10)-5.2)*(-2*(treat=='a') + 3.5*(treat=='b') + 2*(treat=='c')) 
# Simulate binary y to have Prob(y=1) = 1/[1+exp(-L)] 
y <- ifelse(runif(n) < plogis(L), 1, 0) 

##### 1.2 加载数据集至环境-rms 包的要求 ####
# 可以加载 data
data <- data.frame(cholesterol = cholesterol, treat = treat, 
                   num.diseases = num.diseases, age = age, y = y)
ddist <- datadist(data) #
options(datadist="ddist") #  defines data dist. to rms

# 也可以直接加载变量
# ddist <- datadist(cholesterol , treat , num.diseases , age) 
# options(datadist="ddist") #  defines data dist. to rms

##### 1.3 构建含 RCS 的模型 ####
fit <- rms::Glm(y ~ treat + scored(num.diseases) + rcs(age, 4), data = data, x = TRUE)

fit # 模型系数 
# 具有 k 个节点的三次样条将具有 k 个分量，一个常量值（y 截距），
# 一个在建模变量中呈线性的分量（x 值），以及建模变量中的 k-2 个非线性（立方）分量多变的。
extractAIC(fit) 
# AIC 的值，可以改变节点数量，来找到 AIC 最小的节点（即最优模型）
# 但是更多的是要结合临床实际

anova(fit)# 非线性的检验
# p.value_all # 模型的 P value: <.0001
# p.value_nonlinear # 对应的变量的非线性检验: <.0001

# Predict 函数：Obtain predicted values and confidence limits easily 
# varying a subset of predictors and 
# others set at default values

OR.1 <- rms::Predict(fit, age, type="predictions") # 
OR.2 <- rms::Predict(fit, age, type="predictions", fun = exp) # OR 值，Logistics & Cox 模型，需要将 fun = exp，线形回归模型则不需要
OR <- rms::Predict(fit, age, type="predictions", fun = exp, ref.zero = T) # 取中位数为 OR = 1

# ref.zero 参数的解释
# Set to TRUE to subtract a constant from Xβ before plotting so that 
# the reference value of the x-variable yields y=0. 
# This is done before applying function fun. This is especially useful for Cox models to 
# make the hazard ratio be 1.0 at reference values, and the confidence interval have width zero.

# View(OR.1)
# View(OR.2)
# View(OR)
# ddist$limits # 多元回归的复杂性

ggplot(OR.1) # 
ggplot(OR.2) # OR, 经过 exp 的转换后
ggplot(OR) # ref.zero = T，取中位数为 OR = 1


##### 1.4 设置 Reference ####
# 几种策略：1.选择最低点，常见于U型曲线；2）有临床意义的点；3）特殊的分位数；
# OR$age[which.min(OR$yhat)] # 0.6023892

P.10 <- quantile(age, probs = seq(0,1,0.05))['10%'] # -4.79
# P.50 <- quantile(age, probs = seq(0,1,0.05))['50%'] # 6.646342 

ddist$limits["Adjust to", 'age'] <- P.10 # 选择 10 分位数的 age 结果
# ddist$limits["Adjust to", 'age'] <- P.50 # 选择 10 分位数的 age 结果

ddist$limits

# 在新的 Reference 下构建模型
fit <- Glm(y ~ treat + scored(num.diseases) + rcs(age, 4), data = data)
anova(fit)
OR <- Predict(fit, age, fun = exp, ref.zero = T) # 取中位数为 OR = 1
ggplot(OR) # 图形解读，关键在于 OR 为 1 的点

##### 1.5 RCS 图形美化 ####
# 涉及 ggplot 的基本使用，大家可以课后自己优化
# 优化的关键点演示
ggplot() +
  geom_line(data = OR, aes(x = age, y = yhat), linetype = 'solid', size = 1, alpha = 0.7, color = 'red')+
  geom_ribbon(data=OR, aes(x = age, ymin = lower, ymax = upper), alpha = 0.1, fill = 'blue')+
  geom_hline(yintercept=1, linetype=2, color="grey") +
  scale_x_continuous('Age') +
  scale_y_continuous("OR (95% CI)")+
  geom_text(aes(x=15, y=1, label=paste0(
    "P-overall <.0001", # your model pvalue
    "\nP-non-linear <.0001")), hjust=0) + # rcs term pvalue+
  theme_bw() +
  theme(
    axis.line=element_line(),
    panel.grid=element_blank(),
    panel.border=element_blank()
  ) 

##### 1.6 概括：核心 2 步，1）加载数据；2）构建模型, 预测结果 ####
# 1.6.1
dd <- rms::datadist(data)
options(datadist = "dd")

# 1.6.2
# Linear 线性回归样例
res.linear.rcs <- rms::ols(y ~ rcs(x1,4) + x2, data = data, x = TRUE)
pre.linear <- rms::Predict(res.linear.rcs, x, type = "predictions", ref.zero=T)

# Logistics 回归样例
res.logitics.rcs <- rms::Glm(y ~ rcs(x1,4) + x2, data = data, x = TRUE)
pre.logistics <- rms::Predict(res.logitics.rcs, x, fun=exp, type = "predictions", ref.zero=T)

# Cox 回归样例
res.cox.rcs <- rms::cph(Surv(time, status) ~ rcs(x1,4) + x2, data = data, x = TRUE, y = TRUE, tol=1e-25, surv = TRUE)
pre.cox <- rms::Predict(res.cox.rcs, x, fun=exp, type = "predictions", ref.zero=T)



#### 2. 加权线性模型 + RCS (以 Cox 为例) ####
##### 2.1 提取数据模块 ##### 
demo.i <- read_xpt("NHANES_data/2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)
dim(demo.i) # 9971

# 提取生存数据 ：https://www.cdc.gov/nchs/data-linkage/mortality-public.htm

# 连接 National Death Index (NDI) ，补充 NHANES 参与人员的死亡数据
# 得到公开的 linked mortality files (LMF) 提供从调查参与日期到2019年12月31日的死亡率随访数据。
# 公开的数据中仅有随访时间或潜在的死亡原因
# 生存地址原始下载链接：https://ftp.cdc.gov/pub/Health_Statistics/NCHS/datalinkage/linked_mortality/
# 可以下载不同周期的生存数据，下载后进行读取

# setwd("C:/PUBLIC USE DATA")
mort.data_2015_2016 <- NHSMortRead('NHANES_2015_2016_MORT_2019_PUBLIC.dat')
# eligstat
# A value of 1 for ELIGSTAT indicates that the survey participant was eligible for the mortality linkage, 
# a value of 2 indicates the survey participant was under age 18 and not eligible for public release, 
# a value of 3 indicates the survey participant was not linkage eligible due to having insufficient identifying data to conduct data linkage.

# ucod_leading: 10种 潜在死亡原因，001-010，分别对应的原因看课程附件材料
# mortstat: 1生存状态
# permth_int: 从interview到随访时间点，间隔的月份数
# permth_exm: 从MEC-移动检查车到随访时间点，间隔的月份数

dim(mort.data_2015_2016) # 9971
# 转变为大写的列名，和 demo 进行拼接
colnames(mort.data_2015_2016) <- toupper(colnames(mort.data_2015_2016))
# View(mort.data_2015_2016)


##### 2.2 合并数据提取变量  & 去掉 NA 的行 & 衍生变量 ##### 
analyze.sample.data <- demo.i[,c('SEQN',
                                 "RIDAGEYR", "RIAGENDR", "INDFMPIR", "DMDEDUC2",
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]

analyze.mort.data <- mort.data_2015_2016[,c('SEQN', 'ELIGSTAT', 
                                            "MORTSTAT", "PERMTH_INT")]
# 拼接 demo & mort 数据
analyze.sample.data.add.mort <- join_all(list(analyze.sample.data, analyze.mort.data))

dim(analyze.sample.data.add.mort) # 9971    
# View(analyze.sample.data.add.mort) 

# 一键去掉 NA 的行
analyze.sample.data.add.mort.drop.na <- drop_na(analyze.sample.data.add.mort) 
dim(analyze.sample.data.add.mort.drop.na) # 5071

# 衍生低收入 vs 非低收入的变量 PIR.factor
analyze.sample.data.add.mort.drop.na$PIR.factor <- ifelse(analyze.sample.data.add.mort.drop.na$INDFMPIR < 1.3, 1, 0)

# 转换为分类变量
analyze.sample.data.add.mort.drop.na$DMDEDUC2 <- factor(analyze.sample.data.add.mort.drop.na$DMDEDUC2)

##### 2.3 生成复杂抽样 NHANES_design ##### 
NHANES_design <- svydesign(
  data = analyze.sample.data.add.mort.drop.na, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~WTINT2YR,
  survey.lonely.psu = "adjust") # 可以加上 survey.lonely.psu = "adjust" 避免1个PSU报错

summary(NHANES_design)

##### 2.4 复杂抽样的 Cox 回归模型 ##### 
# 死亡和贫困指数之间的相关性
fit <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor + rcs(RIDAGEYR,4) + DMDEDUC2, 
               design = NHANES_design)

##### 2.5 构建含 RCS 的复杂抽样的 Cox 回归模型 ##### 
# 因为 rms 无法直接 svydesign ，因此需要将权重进行转换，变成 rms 包可以处理的形式
data <- fit$survey.design$variables # 从 svydesign 中提取 data
# 添加权重变量至 weights
ori.weight <- 1/(fit$survey.design$prob)
mean.weight <- mean(ori.weight)
data$weights <- ori.weight/mean.weight # 生成用于计算的权重变量，加入到 data 中

dd <- rms::datadist(data)
options(datadist = "dd")

# 得到加权情况下的含 rcs 的 cox 模型
fit.rcs <- rms::cph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor + rcs(RIDAGEYR, 4) + DMDEDUC2, 
                    data = data, weights = weights, normwt = TRUE) 
# if weights are normalization (aka reliability) weights, then normwt should be TRUE

anova(fit.rcs)

##### 2.6 概括：核心 3 步 ####
# 2.6.1 构建复杂抽样下的回归模型
# Cox
fit <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor + rcs(RIDAGEYR,4) + DMDEDUC2, 
                design = NHANES_design)

# 2.6.2 通用-构建含 weight 的数据，并加载至环境
data <- fit$survey.design$variables # 从 svydesign 中提取 data
# 添加权重变量至 weights
ori.weight <- 1/(fit$survey.design$prob)
data$weights <- ori.weight/mean(ori.weight) # 生成用于计算的权重变量，加入到 data 中
# 加载环境
dd <- rms::datadist(data)
options(datadist = "dd")

# 2.6.3 得到加权情况下的含 rcs 的模型(仅是在 1.6.2 的基础上加上 2 个参数，weights = weights, normwt = TRUE)
# 1.6.2
# Linear 线性回归样例
res.linear.rcs <- rms::ols(y ~ rcs(x1,4) + x2, data = data, weights = weights, normwt = TRUE)
pre.linear <- rms::Predict(res.linear.rcs, x, type = "predictions", ref.zero=T)

# Logistics 回归样例
res.logitics.rcs <- rms::Glm(y ~ rcs(x1,4) + x2, data = data, weights = weights, normwt = TRUE) # y 取值为 0，1 
pre.logistics <- rms::Predict(res.logitics.rcs, x, fun=exp, type = "predictions", ref.zero=T)

# Cox 回归样例
res.cox.rcs <- rms::cph(Surv(time, status) ~ rcs(x1,4) + x2, data = data, weights = weights, normwt = TRUE)
pre.cox <- rms::Predict(res.cox.rcs, x, fun=exp, type = "predictions", ref.zero=T)



